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Electromagnetic driving in a honeycomb lattice can induce gaps and topological edge states with a structure 
of increasing complexity as the frequency of the driving lowers. While the high-frequency case is the most 
simple to analyze we focus on the multiple photon processes allowed in the low-frequency regime to unveil 
the hierarchy of Floquet edge states. In the case of low intensities an analytical approach allows us to derive 
effective Hamiltonians and address the topological character of each gap in a constructive manner. At high 
intensities we obtain the net number of edge states, given by the winding number, with a numerical calculation 
of the Chern numbers of each Floquet band. Using these methods, we find a hierarchy that resembles that 
of a Russian nesting doll. This hierarchy classifies the gaps and the associated edge states in different orders 
according to the electron-photon coupling strength. For large driving intensities, we rely on the numerical 
calculation of the winding number, illustrated in a map of topological phase transitions. The hierarchy unveiled 
with the low-energy effective Hamiltonians, along with the map of topological phase transitions discloses the 
complexity of the Floquet band structure in the low-frequency regime. The proposed method for obtaining the 
effective Hamiltonian can be easily adapted to other Dirac Hamiltonians of two-dimensional materials and even 
the surface of a three-dimensional topological insulator. 

PACS numbers: 67.85.Hj; 73.22.Pr; 73.20.At; 72.80.Vp; 78.67.-n 


I. INTRODUCTION 

A topological material or system (e.g., a quantum Hall insu¬ 
lator or a topological insulator) has a bulk gap characterized 
by a topological invariant bearing a non-trivial value ID 12. 
The bulk-boundary correspondence establishes that when in 
contact with the vacuum (or a trivial material) the interface 
between the two media hosts conducting edge states CD- In¬ 
terestingly, the number and chirality of the edge states are 
solely determined by the topological invariants computed for 
the bulk systems. Recently, several studies signaled that topo¬ 
logical edge states can be engineered in an ordinary material 
by applying a time-periodic driving EH3- This sparked the 
interest of diverse communities from graphene ||6l4T2l and re¬ 
lated materials QMI}, to topological insulators GO ED, 
photonic crystals 02, and optical lattices fl8] - [26l , aiming 
to tackle a plethora of issues: characterization of these novel 
edge states 0M2, different signatures in magnetization and 
tunneling (27}? ], the proper invariants entering the bulk¬ 

boundary correspondence (2814301 . their statistical proper¬ 
ties f3~n 132 i. the role of interactions and dissipation (3314351 
and the associated two-terminal (36) J37jj and multiterminal 
(Hall) conductance (35l [388 . 

Floquet theory | [39ll43l is the prevalent tool for the study 
of time-periodic Hamiltonians. Within Floquet theory, the 
solutions of the time-dependent Schrodinger equation can be 
conveniently casted in terms of the solutions of an eigenvalue 
problem in a higher-dimensional space, the so-called Floquet 
space |[39l (431 which is the direct product between the usual 
Hilbert space and the space of time-periodic functions with 
period T = 2tt/CI. The increased dimensionality is at the heart 
of the richness arising in the Floquet quasienergy spectra. No¬ 
tably, when the driving opens a gap between two adjacent Flo¬ 


quet replicas, other replicas (associated to different number of 
photons) develop a hierarchy of ever smaller gaps, each of 
them hosting chiral edge states. The ensuing structure, which 
reminds us of Russian nesting dolls, progressively unfolds as 
higher-order inelastic processes are explored. 

While for high-frequency driving, i.e., of the order of or 
larger than the bandwidth, the system’s stroboscopic evo¬ 
lution (8} can be elegantly described by an effective time- 
independent Hamiltonian EH E2 M, the opposite low- 
frequency regime is trickier to deal with, but might be ex¬ 
perimentally more feasible for many materials like three- 
dimensional topological insulators (431 , graphene mmm, 
or other two-dimensional materials □21. Moreover, it is in 
this regime that the mentioned nesting structure appears and 
the determination of an effective Hamiltonian and the charac¬ 
terization of the associated chiral edge states becomes more 
challenging. 

Here we address the nesting structure of the bulk gaps and 
associated edge states in the Floquet quasienergy spectra of 
honeycomb lattices. To do this we rely on the fact that these 
gaps follow a hierarchy in which the gaps’ widths depend on 
the order of the inelastic processes originating them. This al¬ 
lows us to determine the number of edge states by looking 
first at the largest energy scale (largest gap) and progressively 
moving into the smaller (higher-order) gaps towards the gap 
center. The hierarchy unfolds as new edge states bridge the 
smaller gaps. Honeycomb lattices illuminated by an intense 
circularly polarized laser have attracted much attention in this 
context 03112 but a detailed analysis for frequencies span¬ 
ning both high- and low-frequency regimes is missing. Here 
we provide a systematic derivation of the effective Hamilto¬ 
nians at the crossings between Floquet bands together with 
analytical expressions for the associated contributions to the 
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Chern numbers. 

This work is organized as follows. In Sec. [IT] we present the 
Floquet Hamiltonian for an irradiated honeycomb lattice. In 
Sec.|nl]we discuss the calculation of the Chern numbers of the 
Floquet bands in terms of the low-energy (Dirac) Hamiltonian 
and explain the hierarchy of the corresponding edge states. 
The case of large driving intensity and frequency is analyzed 
in Sec. |IV| where a full map of the Chern number is obtained 
by a direct numerical calculation using the bulk tight-binding 
Hamiltonian. This enables us to show a phase diagram of the 
topological phase transitions for a wide range of frequencies 
and intensities of the driving field. 


II. DRIVEN HONEYCOMB LATTICE 

Let us consider a general system with a Hamiltonian Ho 
(time-independent) in the presence of a time-periodic pertur¬ 
bation V{t). The full Hamiltonian H(t) = Ho + V(t) satisfies 
H(t+T) = H(t), where the period T = 27r/H is determined by 
the driving frequency O. Floquet’s theorem guarantees the ex¬ 
istence of a set of solutions of the time-dependent Schrodinger 
equation of the form \ip a (t)) = exp(-i£ a t/h)\(f> a (t)), where 
| <j> a (t)) has the same time periodicity as the Hamiltonian, 

| <j) a (t + T)) = \4> a (t)) Pf)H43l —this is the equivalent of 
the usual Bloch theorem for systems that are periodic in real 
space. The Floquet states \<f> a (t)) are the solutions of the 
eigenvalue equation HF\<t>a(t)) = e«|</>a(i))- where Hf = 
H - ih !j f is the so-called Floquet Hamiltonian and e a is the 
Floquet quasienergy. 

It is customary, and useful, to introduce the notion of 
the Floquet space, formed by the direct product between 
the Hilbert space and the space of time-periodic functions 
with period T (spanned by the functions e mf2t with n = 
0, ±1,±2,...), so that |</>„(£)} = T,n emUt \ u n)■ When writ¬ 
ten in this basis, the Floquet Hamiltonian Hf is a time- 
independent infinite matrix Hp with copies of Ho in the di¬ 
agonal blocks or Floquet replicas (fixed n). Each diagonal 
block is shifted in energy by nhfl. The time-dependent per¬ 
turbation enters only (if it has zero time-averaged value) in the 
off-diagonal blocks that couple the different Floquet replicas. 

In analogy with the concept of the Brillouin zone for 
Bloch electrons, the quasienergies can be restricted to a 
Floquet zone. Indeed, for every solution 1 4> a (t)) with 
quasienergy e a one can construct another solution 1 4> am (t) } = 
with quasienergy e arn = e a + mhfl, that 
corresponds to the same physical state | ip a (t)). Therefore, the 
eigenvalues are repeated at intervals of hi l and they could be 
restricted to the interval (-ftfI/2, ftfI/2]. While this reduced 
zone scheme is the usual choice, we find it more convenient 
and more insightful, for reasons that will become clear below, 
to work in the extended zone scheme. In that case, to bet¬ 
ter interpret the results, a useful magnitude that complements 
the spectral information (see below) is the time-averaged “lo¬ 
cal” density of states which can be computed as the density of 
states associated to the Floquet Hamiltonian projected on the 


n = 0 Floquet subspace mil] 

Pa(e) = Y i 5 ( £ - £ a) l(aK)| 2 > (1) 

OL 

where \a) is an arbitrary state of the Hilbert space. In the sum, 
the full set of quasienergies e a is kept to ensure that for van¬ 
ishing intensity of the time-periodic potential (and hence of 
the coupling between the Floquet replicas) the original den¬ 
sity of states of the unperturbed system is recovered. Equa¬ 
tion 0 can also be casted in terms of the Floquet-Green func¬ 
tion 1 12.48 1. It is worth noting that recent works point out the 
key role played by the time averaged component of the Flo¬ 
quet eigenstates mw$ ,120, particularly when analyzing 
the transport response of the driven system l38l . 


A. Floquet-Bloch Hamiltonian 

A honeycomb lattice with a single orbital per site can be 
described by the following tight-binding Hamiltonian 

"fttb (t) = Y, e i c \ c i ~ E [ 7 ij(t) c\c 3 + h.c.]. ( 2 ) 

* (W> 

Here c[ and c i are the electronic creation and annihilation op¬ 
erators at site i with energy e,, respectively, and 7 ^ is the 
nearest-neighbors hopping matrix element. We neglect the 
spin degree of freedom throughout this work as it does not 
play any role. 

The effect of the circularly polarized electromagnetic 
field E(t) can be described in a gauge such that E(t ) = 
-(1 /c)dA/dt, where A(t) = A 0 (cos fit, sin fit) is the vec¬ 
tor potential—this describes the situation of normal inci¬ 
dence. Hence, the time-dependent field enters the Hamilto¬ 
nian through the hopping matrix elements (Peierls substitu¬ 
tion): 

t 2n r r t \ 

li'j(t) = Texp^i— A(t) ■ d£J , (3) 

where $0 is the magnetic flux quantum. 

Following a similar procedure as in Refs. j49ll50l we arrive 
at the Floquet-Bloch Hamiltonian, Hf(I^) = E m n H m ,n + 
Sm tn hfll, where TT mjTl = 1 /T/ 0 °° e int ^ n ~ rn ^H(t)dt is the 
(n - to) Fourier component of the time-dependent Hamilto¬ 
nian. Each diagonal block has copies of Ho that account for 
the Floquet replicas; the hoppings between different lattice 
sites within the same replica are the zeroth Fourier compo¬ 
nents of 7 j jit). This is proportional to 7 Jq(z) up to a phase 
that depends on the direction of the hopping, where Jq(x) is 
the zeroth order Bessel function, z = Aoa r 2iv/<l>o denotes the 
field intensity from now on, and a c is the distance between 
nearest neighbors in the honeycomb lattice. This dependence 
on Jq(z) will lead to many interesting behaviors of the topo¬ 
logical characteristics (of any driven lattice) when the inten¬ 
sity reaches a root of Jo(x). The first root at 20,1 - 2.4048 
leads to a topological phase transition that is further explained 
in Sec. [TV] 
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B. Low-energy Hamiltonian 


Close to the Dirac points (K and K' points), the band struc¬ 
ture of the honeycomb lattice is well described by a Dirac 
Hamiltonian, 

T~L{t) = hv F [<j x (k x + j^A 0 cos fit) 

+ s <j y ( k y + ^- c A 0 sin fit)] , 

where up denotes the Fermi velocity, er = ( a x ,a y ) are the 
Pauli matrices for the pseudospin degree of freedom, and s = 
±1 is the valley index. 

For the K valley (s = 1) we obtain the Floquet Hamiltonian 


n%(k) 


/ • 


\ • 


ftfi 


hv F k+ 

0 
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hv F k _ 0 0 

hn ^fA 0 0 

^A 0 C 0 hv F k. 

0 hv F k+ 0 





with k± = k x ± iky. Since the external driving is harmonic, 
and in this approximation it enters linearly in the Hamilto¬ 
nian, only the Floquet replicas differing in ±1 photon will 
be directly coupled with a relative strength ij - ev F A( } /cMl 
[in connection with the lattice Hamiltonian r] = (?>n/2hil)z\. 
Higher-order couplings between two replicas, m and n are in¬ 
direct and of order 0(r} n ~ m \). 


III. HIERARCHY OF DRIVING INDUCED GAPS AND 
EDGE STATES 


The Floquet theory outlined in the previous section enables 
a simple picture of how the driving (in our case circularly po¬ 
larized light) can lead to laser-induced gaps m 00 ED. Here 
we briefly highlight a few points that will be useful later on. 
We start considering the low-energy Hamiltonian of Sec. [11} 
A. For vanishing driving strength, we have the Floquet spectra 
represented in Fig.|T}a) (we take here a projection along a par¬ 
ticular k direction around the K point). The effects of the ex¬ 
ternal driving are expected to be important wherever the Flo¬ 
quet replicas corresponding to different values of the Fourier 
index n become degenerate. This happens at half-integer mul¬ 
tiples of hill2. In Fig. [T}a) the crossings at £o = 0 and 
e 1 !2 - hil/2 are marked with gray circles. Interestingly, for 
circularly polarized light all these degeneracies are lifted (in¬ 
cluding the degeneracy between the bands with n = 0 at e = 0) 
with different strengths. In the low intensity limit (// « 1), the 
magnitude of each anticrossing [of order 0(r] An )] is ruled by 
the difference An among the associated replicas, thereby es¬ 
tablishing a hierarchy. This is schematically represented in 
Fig. [3(b) and (c). 

Once the degeneracies develop into gaps, something inter¬ 
esting in the physics of topological systems happens: Edge 
states develop within each anticrossing and these states can 
co-exist with the continuum spectrum provided by other Flo¬ 
quet bands (these bands also have a gap of smaller width). 
The chirality and the robustness to disorder of such Floquet 



FIG. 1. (a) Sketch of the dispersion of the first replicas around n = 0. 
The crossings occur at the Floquet zone center, e = 0, and at the 
Floquet zone borders e = ±hfl/2 are depicted with circles. Note 
that for e = 0 the crossings involve replicas where m + n = 0 and 
are of order rj n ~ rn ' ([ n - m\ even), while in the case e = hkl/2 the 
crossings involve replicas where m + n = ±1 and are of the order 
(| n _ m | odd), (b) and (c) Cartoon representations of different 
crossings for e = 0 and hfl/2 respectively, ordered hierarchically 
according to their magnitude. The special case of the first doll in (b) 
represents the anticrossings at the Dirac point of the n = 0 replica. 
This occurs because of a second-order process involving the emission 
and reabsorption of a photon. 


edge states were explicitly shown in Ref. tED and more re¬ 
cently, other authors pointed out that this could be a general 
fact also in time-independent systems |[52l . In the following 
we will exploit the structure shown in Fig. [T]to systematically 
and progressively unfold our Floquet Russian nesting doll. At 
each step we will obtain an effective Hamiltonian describing 
the corresponding anticrossing, and the number and chirality 
of the edge states bridging it. The latter requires the deter¬ 
mination of the relevant topological invariants that we briefly 
discuss in the next subsection. We then follow with our results 
for the low- and high-frequency regimes. 


A. Topological invariants for Floquet bands 

The Chern number associated to a given Floquet-Bloch 
band a is given by 

C a — — (j) {Uak\^k\^ak) ' dk 

l C r , (6) 

= Im / {^k v 'Uctk \&k x 'U'Oik ) d 
7r J BZ 
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where \u a k) is the periodic part of the Bloch eigenfunction 
and C is the contour of the Brillouin zone (BZ) j53j. Alterna¬ 
tively, Eq. 0 can be cast in the form 

C a = ±[ T ak - dS k , (7) 

Z7T JBZ 

with 

T V—' ( Wafc I Vfc-ff/c | Vfc-fffc ) 

r afc = Im 2, - 7 -tp-, (8) 

/3*a ( £ afc e /3fc)“ 

where I’„fc is the Berry curvature. The peaks in the Berry 
curvature that occur at the points in the BZ where the bands 
are quasi degenerate yield the main contribution to C a . If the 
curvature decays fast enough (which happens when rj -* 0) 
the sum of these contributions is the exact calculation of C a . 
We will make use of this fact in Sec. IIIIBI where an effec¬ 
tive Hamiltonian is derived for the quasi degenerate subspace. 
We also note that though the topological invariants may seem 
very abstract objects they have recently been measured in cold 
matter experiments E3- 

In periodically driven systems to accurately account for the 
edge states one must rely on the winding number W (s) due to 
the infinite periodicity of the Floquet spectrum Il4l f28l . When 
the winding number is evaluated, inside a gap counts for the 
net number of chiral edge states. In connection with the Chern 
number, the difference of winding numbers evaluated at ener¬ 
gies enclosing a band, yields the Chern number of that single 
band. From now on, we will only need to evaluate the wind¬ 
ing number in the two distinct Floquet band gaps, the gap at 
the center of the Floquet zone [W(e o) = W(0)], and the gap 
at the edge of the Floquet zone [kF(ei/ 2 ) = W(hti/ 2)]. This 
topological invariant can be obtained in terms of the evolution 
operator but, here we use an alternative approach proposed 
in Ref. l28l . that consists in truncating the Floquet Hamilto¬ 
nian between the replicas —M and M up to a sufficiently large 
M (note that each extra replica adds two bands to the Floquet 
spectrum). The difference between the number of chiral edges 
states between the a and the (a + 1) Floquet bands will be 
given by 

W(e a ) = £ Cp, (9) 

/3=-(2M+l) 

for a quasienergy e a inside the gap, provided that enough 
Floquet replicas are counted until the sum converges. This 
happens when taking a larger M leaves W'(e a ) unchanged, 
meaning that all relevant crossings between different replicas 
are included in the Floquet zone. Notice then that the con¬ 
tinuum Dirac model is only appropriated as an approximation 
and requires a finite number of replicas. 

A direct evaluation of Eqs. ([7J and ([8]) usually requires the 
use of numerics and the highly peaked Berry curvature ren¬ 
ders the calculation easier for high frequencies. In this regime 
we can characterize the topological properties of the Floquet 
bands and the corresponding edge states using the bulk Flo¬ 
quet Hamiltonian, as seen in Sec.|IV| 


B. Multiple photon processes for low-frequency driving 

In this section we will apply a consistent method to obtain 
the number of edge states inside the driving induced gaps for 
the particular case of the honeycomb lattice. To do this we will 
take advantage of the hierarchy of these gaps, which scale as 
a power of rj with the exponent being the number of photon 
processes. 

To obtain the winding numbers W(eq) and W(e 1 / 2 ) asso¬ 
ciated to the driving induced gaps at the Floquet zone center 
and at the Floquet zone edge, respectively, we must calcu¬ 
late the Chern numbers of all the Floquet bands below them. 
As outlined by Eq. ([8]) the main contributions to the Chern 
number of each band comes from the points in the k space 
where the energies are nearly degenerate. For a vanishing in¬ 
tensity the degeneracies will appear at the crossings of the Flo¬ 
quet replicas. When turning the electromagnetic field on, all 
the degeneracies will be lifted, opening gaps at every avoided 
crossing. 

Fet us use the limit of vanishing intensity to calculate the 
Chern number C a of the a band. This can be obtained as the 
sum of all the contributions from the fc-space regions where 
an avoided crossing occurs. We will denote the contribution 
coming from a point k p c , where the a band has avoided cross¬ 
ings with the (a + 1) band as c^ p a , and if an avoided crossing 
occurs at a (possibly different) point k p : a with the (a - 1) 
band it will be denoted by c ] °, w a . So the sum that yields the 

Chern number is C a = Y. p c p p a + Z,/ C P V*' Here eac ^ con " 
tribution is obtained from Eq. 0 integrating only near the 
avoided crossing -for any finite intensity this is an approx¬ 
imate result but taking the limit where the intensity goes to 
zero the calculation becomes exact. 

Since each avoided crossing means a contribution to the 
Chern number for the bands above and below it with opposite 
signs (c£ft a = ~Cp P a ), when adding all the Chern numbers up 
to the a band to obtain the winding number, most of these 
contributions will cancel out except for the lasts ones (note 
that the first band in the truncated Floquet spectrum has no 
band below and no crossings c^ 0 ™; see Fig. pi). So we obtain 
W a = Z p e/ 1 /.. Since these contributions are the only ones 
that determine the number and chirality of the edge states we 
can drop the superscript in the following. 

We can see in Fig. |T]that the degeneracies appear at k p{) - 
2pko for the gap at the Floquet zone center (eo), and k p , 1/2 = 
(2 p + l)jfc 0 for the gap at the Floquet zone edge (£ 1 / 2 )- being 
t’P ko = f>/2 and p being an integer number. In order to get 
the contribution near the anti-crossing between two replicas 
it is sufficient to derive a 2 x 2 effective Hamiltonian, valid 
close to k p p, with B either zero or one-half. By writing this 
Hamiltonian as 

nf(k,p,/3) = v F h Pj p(k)-(T + £pI 1 (10) 

one can obtain the contribution to the Chern number by calcu¬ 
lating 

Cpp — ^ * (ilp, x dkykip.p) d k. 


( 11 ) 
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with h P: p = h Pi p/\h Pt p\. 

To obtain an explicit form for 'H C F (k,p, ft) we start by 
making a unitary transformation of the pseudospin basis. The 

basis {l/\/2, ±exp(*0fc)/\/2| diagonalizes every diagonal 
block in Eq. ([5} (Floquet replica) describing the Dirac cone 
with eigenvalues shifted by nhQ for the n-th Floquet 

replica. As depicted in Fig.|T]ieplicas indexed by m and n will 
cross at £ 0 when m + n - 0, while the crossing will occur at 
e 1 / 2 if to + n = 1. Hence, we must calculate the effective cou¬ 
pling between the replicas (-m) and (m) or (to + 1) according 
to whether we are evaluating c m p or c m p / 2 . This is achieved 
by a standard procedure based on the projected Green’s func¬ 
tion (or decimation procedure). Namely, if Gj?(u, k) denotes 
the Floquet Green’s function, Gf(w, k) = [col - Tip (fc)]^ 1 , 
we define the effective Hamiltonian, in this case, as 

nf(k, P , p) = n° F (k, P , p) - G^iphfi, k p0 ), (i2) 


where G F (co, k Pt p) = P^ pG F (co,k Pt p)P Pt p is the Green’s 
function projected on the degenerate subspace (n + m = 0 
or m + n - 1) and evaluated at the crossing k Ptl g, P p4 3 is 
the corresponding projector operator, and 'H a F {k 1 p, /?) is the 
projected Floquet Hamiltonian in the absence of the radiation 
field. Excluding the special case of the crossing at fco.o treated 
later, one readily finds, to the lowest non-trivial order in 77 (see 
note li55l ), that 

h p p(k) = r] s va p ko [cos (s p 6k) x + sin (s p 9k) y] 

+ (-k + s p ko)z (13) 

where s p = 2 p, s^ 2 = 2p + 1, tan0fc = k y /k x and a& is a 
numerical factor (see the Appendix). In order to calculate the 
contribution to the Chern number, Eq. it is necessary 
to transform back to a fc-independent basis since the unitary 
transformation we used depends on 9^. This implies a rotation 
of the effective field h p p = Il( 0 k ) h p p. Using polar coordi¬ 
nates we have. 


1 r _ _ _ p 

Cp/3 = — / h p @ • (dghpfi x G/.; hp.£ ) ~ dk dO 

47T J \h fl J 


1 

47T 


1 

\KA Z 

hp,(3 * (dohpfi x dkhp,p') 1 ^ dk dO 


1 r 1 

— / hp^-tR^dgRhp'pxdkhp^)— - -dkd 6 . 

4:77 J 

(14) 


Notice that we took advantage of the fast convergence of the 
integrands and extended the integration to the entire k space. 
The last integral gives zero for h Ptl 3 of the form of Eq. {Il- 
while the other can be done explicitly to obtain 


/ 



Retaining the lowest order in 77 consistent with the approxi¬ 
mation made to obtain h Pjl g, we get 


C-pp - 


(16) 


This is one of our central results. The same derivation can be 
obtained for the expansion around the K' valley, and the total 
contribution (up to the proper order) to G a is twice one 
per each valley. The equality ( fl6| ) could have been anticipated 
from Eq. ( fl3| if one recalls that c p j 3 is related to the number of 
times hpj--i(k) winds around the Bloch sphere as k explores 
the Brillouin zone. The angular dependence of h p p(k) is 
related to the effective coupling between the two degenerate 
replicas through the intermediate ones. From the decimation 
procedure one can infer that the factor in the angular depen¬ 
dence equals the number of replicas decimated plus one or, in 
other words, it is the difference between the Floquet indices of 
the two replicas involved in the avoided crossing. The latter 
makes clear that |s^| is the order of the photon processes that 
lead to the avoided crossing . Following this algorithm when 
looking at the next crossing, p + 1 , the involved replicas will 
be +2 replicas apart, so s ^ +1 = + 2 . 

The only exception to this rule is the particular case of cq.o 
which only comes from the renormalization of the to = 0 
replica and there are no intermediate replicas involved. In this 
case we have 


^o.o(fe) = -2p 2 k 0 x + kz. (17) 

It is clear from the above expression that the value of co.o is 
determined by the last integral in Eq. • 0 - leading to 

co,o = - ^ ■ (18) 

Since we must count both Dirac cones (K and K' valleys) to 
get the total contribution to the Chern number, we get a total 
of -1 for the edge state connecting the K and K' valleys. 
This is the only case where a contribution with a minus sign 
is observed and interestingly enough is a contribution where 
the process involved is of the same order of c-| 0 = 2. This 
allows the edge states of the two K and K' valleys to mix 
with each other and makes a total of 2(co,o + ci^) = 3, which 
is compatible with what is observed in Figs.[2](a) and (c). 

Figure [2] depicts the averaged local density of states near 
the edge of a semi-infinite plane for the radiated honeycomb 
lattice, using the recursive Green’s-function method described 
in Ref. Q3. Here we can also observe the higher-order gaps. 
Since the width of the gap is of order we use a loga¬ 

rithmic scale expanded around £0 = 0 in Fig.[2](c) and around 
£172 = hfl/2 in Fig. [21(d). This allows us to zoom in the spec¬ 
trum up to a cut-off quasienergy denoted by e. This threshold 
is imposed arbitrarily, but constrained by the number of con¬ 
sidered replicas and numerical precision. Note also that the 
weights of different replicas decay exponentially as 77” 1 for 
the TO-th replica; this is evident from the logarithmic scale in 
the color bar of Fig. [2] 

The procedure presented in this section accounts for the 
firsts orders of the generation of gaps and edge states and also 
has the advantage of retaining the largest gaps and the pri¬ 
mary contributions to the averaged density of states. This pro¬ 
cedure is correct if the quasienergies of the replicas involved 
lie within the van Hove singularities of each replica; other¬ 
wise, deviations due to the inaccuracy of the low-energy Dirac 
Hamiltonian appear and a full tight-binding model is required. 
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FIG. 2. (Color online) fc-resolved local density of states near the edge of a semi-infinite honeycomb-lattice plane (in logarithmic scale in 
color). The plane is irradiated with hQ = O.27 and z = 0.05. In (a) the gap at eq is shown for both the K and K' valleys; in (b) the gap at e \/2 
is shown near the K valley only. In the lower panel (c) the energy scale is expanded exponentially around so = O7 up to a minimum cutoff 
energy e ~ 3 x 1CT 8 7. In panel (d) the energy is expanded exponentially around ej/ 2 = O.I7 up to a minimum cutoff energy E - 6 x 1CT 10 7 
[meaning that the interval (-£, i) is not shown]. The lower panels (c) and (d) show the nested hierarchy in powers of 7 of the developed gaps 
and their edge states. 


While there is a plethora of edge states appearing inside the 
gaps, some states might not be measurable simultaneously. In 
a transport experiment with non-irradiated leads only those 
which contribute significantly to the time averaged density of 
states will give a transport channel at the edge of the sam¬ 
ple. In the approach of small r] the main contribution to the 
time-averaged density of states will be given only by the first- 
order gap and its associated edge state at e ~ Ml/2, and to the 
second-order gap for e ~ 0. For more details on the conduc¬ 
tance for a transport calculation we refer the reader to f38l . 


C. High intensity and high-frequency driving 

Now, let us briefly comment on the regime of high frequen¬ 
cies. Because of the reduced number of inelastic processes 
imposed by the higher energy cost, this regime is naturally 
less complex than the one addressed in the previous section. 
Notwithstanding, other difficulties must be taken care off. In¬ 
deed, for frequencies comparable to the band width, the low- 


energy approximation does not hold and the full tight-binding 
Hamiltonian is better suited in this case. For low intensities 
the system can still be solved perturbatively in the Floquet 
space or exactly for the truncated Floquet Hamiltonian, taking 
care of including at least all the Floquet replicas that fit in the 
replicas bandwidth, namely, A/hCl, where A is the bandwidth 
[in our case A shrinks as 6 ^/Jq(z), where z is the driving in¬ 
tensity, see Sec II A|. 


As the driving intensity is increased, higher-order inelastic 
processes are reinforced. Consequently the solutions for the 
infinite Floquet Hamiltonian are spread among more Floquet 
replicas. To obtain a numerical solution we truncate the Flo¬ 
quet Hamiltonian between the -M and M replicas. We must 
include as many replicas as needed for the winding number to 
converge. For example in Fig. [3] even though A/hfl < 3, we 
need five Floquet replicas to obtain the correct result. 


The construction of the winding number is also depicted 
in Fig. [3] where each Floquet band has its associated Chern 
number at the left side, and the two relevant gaps at e = 0 
and htt/2 have their associated winding numbers. The en- 
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FIG. 3. Density of states near the edge of a semi-infinite honey¬ 
comb lattice. The lattice is driven under an electromagnetic field 
parametrized with frequency hfl = 1.87 and intensity z = 1.2. 


hancement of the inelastic processes may lead to unexpected 


topological phase transitions as discussed in Sec. IV 


IV. TOPOLOGICAL PHASE TRANSITIONS 


number to converge one must look at the replicas that would 
reach the £ () and e \/2 points for vanishing intensity —for 
higher intensities more replicas are needed as explained be¬ 
low. 

The bandwidth of the ?i-th replica lies between nhQ± 37, so 
it will be an overlap of different Floquet bands at frequencies 
hCl < 37 /n for the £0 crossing , and at hfl < 3"f/(n - for 
the e 1 /2 crossing (assuming n > 1 and vanishing intensity). 
This behavior is shown in Fig. [4] where for low intensities 
a topological phase transition occurs every time a new pair 
of replicas enters in the description of the system. For low 
intensities, the bandwidth of the replica shrinks proportionally 
to 7Jo(z) ~ 7(1 - z 2 ), which can be seen as down going 
parabolas at htt = 37 /n in Fig. Wha) and at ftfl < 3"f/(n - |) 
in Fig. 0(b). 

The above deduction is based on the fact that for low inten¬ 
sities, the hoppings between one site in the n-th replica and 
one site in the m-th replica are proportional to 7 J n - m (z) ~ 
7 z n ~m . This means that for low intensities the dominant cou¬ 
pling is the zeroth-order one, i.e., the one within the same 
photon subspace. As one increases the intensity this assump¬ 
tion no longer holds and the coupling between neighboring 
replicas can achieve larger values, forcing the eigenfunctions 
that solve the Floquet Hamiltonian to be spread among many 
Floquet subspaces (replicas). For higher intensities the effects 
of introducing a new replica in the calculation extends beyond 
the replica’s bandwidth and to correctly address the topology 
of the system one must include a larger number of replicas in 
the calculation. This is the explanation for the lines with a 
positive slope in Fig. [4] that mark a topological phase transi¬ 
tion. 


In the previous section we showed that for low frequencies 
there is a growing number of edge states as a larger number 
of replicas are included in the calculation. There is, how¬ 
ever, a natural limitation to this procedure, when the k • p 
approach no longer describes correctly the topology of the 
Floquet bands involved. In the case of the honeycomb lat¬ 
tice the van Hove singularity sets this energy threshold. The 
van Hove singularities lie at energies of nhtt ± 7 for the n-th 
Floquet replica so it will be well described at the £0 crossing 
only for frequencies such that hfl < 7/n, and at the £-72 cross¬ 
ing for hfl <7/(n - |), assuming n > 1. To illustrate this let 
us choose a frequency of Ml = O.O57. In this case the replicas 
from to = -20 to 20 are well described at £0, and the replicas 
to = -18 to 19 are well described at Ey 2 . So for low fre¬ 
quencies and moderate amplitudes the low-energy approach 
ensures to take into account all relevant Floquet replicas nec¬ 
essary for the calculation to converge, and to accurately ad¬ 
dress the number and chirality of the edge states relevant for 
transport. 

As we increase the frequency, this number rapidly drops 
and one must use the full tight-binding Hamiltonian to de¬ 
scribe the bands in a wider energy range. To find out how 
many replicas are needed for the calculation of the Chern 


Another interesting behavior is the transition at z = 2.4048 
for all frequencies, marked by a vertical line in Fig. [4](a) [and 
less resolved in (b)]. At this point the hopping between sites 
that belong to the same replica vanishes; this is the first zero 
207 of the Bessel function Jo(x). Some numerical noise can 
be seen in panels (a) and (b) for low intensity because of the 
vanishing width of the highest-order gap considered, and also 
there is noise at some lines depicting a topological phase tran¬ 
sition since the gap closes at every phase transition. The cal¬ 
culation time rapidly grows as more replicas are considered, 
which is the reason for the blank slices in the bottom left of 
(a) and (b). For larger values of z a quasi periodic pattern 
is observed due to the Bessel’s functions quasi-periodicity. 
This regime is not shown here because the intensities involved 
are extremely high for a possible experimental realization and 
for the assumptions made when modeling the electromag¬ 
netic field, and the system can become unstable against slight 
changes from circularly to elliptically polarized light, as stud¬ 
ied in <471 for high frequencies. Instead of the winding num¬ 
ber, a map of the Chern number is presented in 137ft . Besides 
that, some phase transition could remain hidden for the time 
averaged transport in a multiterminal scattering configuration 
138ft . since the corresponding edge states could bear no weight 
in the time-averaged density of states. 
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FIG. 4. Map of the winding number W (e), calculated with the full Floquet-Bloch bulk Hamiltonian, for eq = 0 in panel (a) and £ 1/2 = hfl/2 
in panel (b). A maximum number of 11 replicas has been used throughout (M = 5). 


V. CONCLUSIONS 

Characterizing the topological properties of driven sys¬ 
tems in general, and honeycomb lattices in particular, is 
crucial for many studies pursuing novel Floquet topological 
phases EE EH ED. In this paper we address the calcula¬ 
tion of the topological invariants in a wide range of parame¬ 
ters, from high to low frequencies. The Floquet quasienergy 
structure becomes progressively more complex when the fre¬ 
quency becomes much smaller than the bandwidth. In par¬ 
ticular, within a small photon-energy range we find a nested 
structure of gaps of different widths, which are proportional to 
a power of the electron-photon coupling, the exponent being 
related to the order of the inelastic processes. 

Interestingly, Floquet edge states develop within each gap 
even in the presence of a continuum of other Floquet bands 
provided that the edge states and the continuum have very dif¬ 
ferent spectral weights among the replica’s subspace. This 
allows one to devise a scheme for the determination of the 
number and chirality of the edge states where this information 
is progressively obtained as higher-order inelastic processes 
are included. This procedure is limited by the ratio between 
the system’s bandwidth and the driving frequency. 

The first stage of the scheme presented here is the calcula¬ 
tion of an effective Hamiltonian, which is done analytically. 
This effective Hamiltonian is aimed at describing the Floquet 
quasienergy structure, rather than the time evolution, and al¬ 
lows one to compute the topological invariants in a broad set 
of driving frequencies and intensities. For low frequencies 
we have derived the contributions to the Chern numbers, con¬ 
structively matching the numerical results obtained using re¬ 
cursive Green’s functions for the Floquet-Bloch tight-binding 
Hamiltonian. 

For higher frequencies and a vast set of intensities the nu¬ 
merical evaluation of the winding number is summarized in a 
map of topological phase transitions. The main features are 
the lines that mark a topological phase transition where dif¬ 
ferent numbers of Floquet replicas become degenerate. This 


allows one to tune the radiation parameters in order to obtain 
a specific number of edge states. 
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Appendix: Effective Hamiltonian 

In this section we deduce the effective Hamiltonian that de¬ 
scribes the crossing of different Floquet replicas, say the repli¬ 
cas labeled by m and n. The crossings at £0 occur for to = -n 
and the crossings at £^2 occur for m = —n + 1, where n > 1. 
The most simple way to evaluate the effective coupling of the 
replicas is to make a change of basis that diagonalizes each 
subspace of the Floquet Hamiltonian 'Hp(k) in Eq. ([5]) to get 


••• 



- H^(k) 

V(0 k ) 

0 

... v(e k y 

H ( 0 °\k ) 

V(9 k ) - 

0 

v{e k y 

H^ik) ••• 

: 


: •• 


(A.l) 
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where 


HQ n \k) = u F ^g _° fc j + nv F 2k 0 I 
V(O k ) = r,k 0 e~^ ( _^_ iX j , 


(A.2) 


where e lA is the trivial phase between the two basis vectors 
(in the following valued in A = 0), 9 k is the angle between k 
and the x-axis, and 2 v F ko has replaced hi 1 to make evident 
the crossing point in the k space. As stated in Sec. IIIB the 
values of the modulus of k where the replicas cross will be 
k p j 3 = 2ko(p + [3 ), where j3 - 0 or 1/2, according to whether 
we are looking at the crossings at Sq or £ 1 / 2 , respectively. The 
next step is to apply the decimation procedure, thus eliminat¬ 
ing the replicas in between, to renormalize the effective hop¬ 
pings that couple the desired replicas (the renormalization of 
the diagonal terms is irrelevant, for the purposes of calculating 
the Chern numbers, and will be neglected for simplicity). In 
the case that the replicas to = 0 and n = 1 it is straightforward 
to see from Eq. (|A. 1 [> that the effective Hamiltonian will be 


Hf (£, 0 , 1 / 2 ) = w F 


-k + 2ko -rjk oe l6k \ 

- 77 k 0 e zek k J ’ 


(A.3) 


since there are no replicas in between to decimate. At the 
same energy e 1 / 2 the next crossing will occur for the replicas 
to = -1 and n - 2 , and we will need two steps of decimation 
for the replicas zero and one. 

Then the decimation of two replicas will have the effect of 
accumulating two orders more in the coupling strength and 
in the phase factor, resulting in a coupling proportional to 
Tj 3 e~ 3iek . The calculation can be performed to easily obtain 
the effective Hamiltonians. Expressed in terms of h p p(k) 
the calculation yields, 


'H c p (k,p,(3) = v F h p p(k) ■ a + EpI. Using this expressions 
we can evaluate Eq. ( fl4| ) to calculate the contribution of these 
crossings to the winding number, i.e., c 0 i i /2 = 1, C \AI 2 = 3, 
c 2 j i / 2 = 5, etc. 

The same procedure can be applied to the crossings at £q, 
starting from the crossing of the replicas m = -1 and n - 1 
at kifi, where only one decimation step is needed, giving a 
effective coupling proportional to n 2 e -2 “ fc . The next crossing 
of the replicas m = -2 and n = 2 will accumulate two orders 
more in these factors, and so on. The explicit calculation gives 


hip(k) = if£ 2 k 


h 2 p{k) = r / 4 £4 £ 

^3,o(fc) = if L k 


1,0 2 + ( _ £ + ^ 1 , 0 ) z 

2,0 — + {-k + k’ 2 , 0 ) Z 
JLo 


3,0 3200 + + ^ 3,o) ^ 


(A.5) 


where h p> p(k) is defined as before. It is straight forward to 
see that c\p = 2, c 2 p = 4, = 6 , etc. 

The only exception is the calculation of hop, which has 
been already addressed by Oka and Aoki (3). This time, the 
effective Hamiltonian is the renormalized Hamiltonian of the 
to = n = 0 replica which has a crossing of its own bands at 
the Dirac point hop = 0. The degeneracy is lifted due to the 
coupling with the replicas ±1, and the effective Hamiltonian 
is described by Eq. 0 . and can be equally expressed as, 

hop(k) = -rf 2koX + kz . (A. 6 ) 


In this case, it is important to rotate back to a /c-independent 
basis as explained in the text. This is done with the rotation 
matrix 


ho,i/ 2 (k) - ~ vii ko,i /2 + {~k + £ 0 , 1 / 2 ) ^ 

3 

hi tl / 2 (k)= -rf £3£1,1/2 — + (-£ + £1,1/2) z 

80 

^ 2 ,l/ 2 (£) = -rf £5 £2,1/2 7^ + (-£ + £2,1/2) Z 

(A.4) 

where the unit vector = cos ( n0 k ) x + sin ( nd k ) y winds 
n times in the xy plane around the z axis as we move £, and 


/ 0 sin 0 cost? ' 
R(9) = I 0 -cos 9 sind 

\i 0 0 

that satisfies the following useful identity 


R 1 dgRh p p = x x h p p 


(A.7) 
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